//
//  Simulation.cpp
//  phy_four
//
//  Created by Snir Gazit on 3/12/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Simulation.h"

typedef boost::chrono::system_clock::time_point TimePoint;
typedef boost::chrono::duration<double, boost::ratio<3600>> Hour;
typedef boost::chrono::duration<double, boost::micro> MicroSec;

double calcWR(double mu, double U, int m) {
	double c1, c2, c3, md = double(m), arg;
	arg = 0.25 * pow(mu, 2) / U;
	c1 = 0.25 * pow(U, -(3. / 2.) - double(md) / 2);
	c2 = mu * gsl_sf_gamma((3. / 2.) + md / 2)
			* gsl_sf_hyperg_1F1((3. / 2.) + md / 2, 3. / 2., arg);
	c3 = sqrt(U) * gsl_sf_gamma(1 + md / 2)
			* gsl_sf_hyperg_1F1(1 + md / 2, 1. / 2., arg);
	return c1 * (c2 + c3);
}
void calcWeight(double muIn, double UIn, double g, int maxBond,
		ostream &outFile, BoltzWeightType &Rweight) {
	double mu = muIn * g, U = UIn * g;
	double n, c, w;
	c = calcWR(mu, U, 0);
	for (int i = 1; i < MAX_BOND; i++) {
		n = calcWR(mu, U, i);
		w = (n / c);
		c = n;
		Rweight[i - 1] = w;
	}
}
// WormLoop Algorithm

// Worm Alg
void OnModelLoop::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
// OutFile
	string runNum = qmcParameters.runNum.getVal();
	// OutFile
	String fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Starting N=4 Simulation - Run Number - " + runNum << endl;
	// Calc Weights
	qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzRatio_.fill(1.);
	else
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(),
				qmcParameters_.g.getVal(), qmcParameters_.maximalBonds.getVal(),
				outFile_, boltzRatio_);
	// Init Lattice
	boost::array<int, DIM> sizes;
	sizes.fill(qmcParameters_.L.getVal());
	if (qmcParameters_.simType.getVal() == "HIGGS")
		sizes[DIM - 1] = qmcParameters_.T.getVal();
	lat_.init(sizes);
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	CreateLoopMoveXY_.init(&eng_, &lat_, qmcParameters_, boltzRatio_);
	chooseLoopType_.init(&eng_, BinomialVariableParams(1, 0.5));
// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzRatio_, &outFile_);
// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	measures_.updateMinNumOfMeasure(totalNumberOfMeasure);
}
void OnModelLoop::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	for (int i = 0; i < 10 * numOfThermal; i++)
		nextMove();
	double cLengthXY = ba::mean(wormLengthXY_), nLengthXY;
	bool conv;
	do {
		wormLengthXY_ = acc();
		for (int i = 0; i < numOfThermal; i++)
			nextMove();
		nLengthXY = ba::mean(wormLengthXY_);
		outFile_ << "Worm Length XY - " << nLengthXY << endl;
		double ratio = fabs(1. - (cLengthXY / nLengthXY));
		conv = ratio < 0.01;
		cLengthXY = nLengthXY;
	} while (!(conv));
//CreateLoopMove_.setAvWormLength(2 * ba::mean(wormLength));
	outFile_ << "Final Worm Length XY - " << ba::mean(wormLengthXY_) << endl;
	double wormLength = ba::mean(wormLengthXY_);
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / wormLength);
//      closedRate = 1;
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

OnModelLoop::~OnModelLoop() {
	outFile_.close();
}

void OnModelLoop::simulatedAnnealing() {
	outFile_ << "Starting Simulated Annealing" << endl;
	double gHighT = 1. / qmcParameters_.gHighT.getVal();
	int numOfSMA = qmcParameters_.numOfGSA.getVal();
	double dg = (qmcParameters_.g.getVal() - gHighT) / double(numOfSMA - 1);
	double nLengthIsing, nLengthXY, nTotal;
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	vector<double> gSA;
	for (int i = 0; i < qmcParameters_.numOfGSA.getVal(); i++)
		gSA.push_back(gHighT + dg * double(i));
	for (auto it_g=gSA.begin();it_g!=gSA.end();it_g++) {
		double cg=*it_g;
		setCoupeling(cg);
		outFile_ << "Current g -" << cg << endl;
		wormLengthXY_ = acc();
		for (int i = 0; i < numOfThermal; i++)
			nextMove();
		nLengthXY = ba::mean(wormLengthXY_);
		numOfThermal = int(
				120 * double(lat_.getNumOfSites()) / double(nLengthXY));
		outFile_ << "Worm Length XY MA - " << nLengthXY << " numOfThermal "
				<< numOfThermal << endl;
	}
	numOfThermal = int(5000 * double(lat_.getNumOfSites()) / double(nLengthXY));
	wormLengthXY_ = acc();
	for (int i = 0; i < numOfThermal; i++)
		nextMove();
	nLengthXY = ba::mean(wormLengthXY_);
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / double(nLengthXY));
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}
void OnModelLoop::setCoupeling(double g) {
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzRatio_.fill(1.);
	else
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(), g,
				qmcParameters_.maximalBonds.getVal(), outFile_, boltzRatio_);
	CreateLoopMoveXY_.setCoupeling(boltzRatio_, g);
}

bool OnModelLoop::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	wormLengthXY_ = acc();
	while (measures_.checkNotDone()) {
		nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump();
		}
		if (endTime) {
			this->dump();
			break;
		}
	}
	outFile_ << "Final Worm Length XY - " << ba::mean(wormLengthXY_) << endl;
	outFile_ << " Erase Prob - " << CreateLoopMoveXY_.getEraseProb() << endl;
	return !endTime;
}
void OnModelLoop::nextMove() {
	CreateLoopMoveXY_.setLoopType(chooseLoopType_.draw());
	wormLengthXY_(CreateLoopMoveXY_.makeMove());
}
void OnModelLoop::dump() {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath);
	lat_.dump(outPath);
	int totalBonds;
	eng_.dump(outPath);
	String fileName = outPath + "/sim_params.xml";
	// Write to XML
	pugi::xml_document doc;
	pugi::xml_node rootNode = doc.append_child("SIM_PARAMS");
	// Parameters
	pugi::xml_node node = rootNode.append_child("CLOSED_RATE");
	std::string data = boost::lexical_cast < std::string
			> (measures_.getClosedRate());
	node.append_child(pugi::node_pcdata).set_value(data.c_str());
	doc.save_file(fileName.c_str());
}

void OnModelLoop::read() {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	outFile_ << "Measures Loaded"<< endl;
	lat_.read(outPath);
	try {
		int totalBonds;
		bool testConf = lat_.testConfiguration(totalBonds);
		outFile_ << "Total Bonds Calculated - " << totalBonds
				<< " total Bonds previous sim - " << lat_.getTotalBonds()
				<< endl;
		if (totalBonds != lat_.getTotalBonds())
			throw "Total Bonds after Read don't match";
		if (!testConf)
			throw "Bad Load of Configuration";
	} catch (string &e) {
		outFile_ << e;
		exit(1);
	}
	outFile_ << "Lattice Loaded"<< endl;

	eng_.read(outPath);
	String fileName = outPath + "/sim_params.xml";
	// Read From XML
	pugi::xml_document doc;
	try {
		pugi::xml_parse_result result = doc.load_file(fileName.c_str());
		if (!result)
			throw result;
	} catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	pugi::xml_node rootNode = doc.child("SIM_PARAMS");
	measures_.updateClosedRate(readXmlData<int>(rootNode, "CLOSED_RATE"));
}
